# engine mean value model only calculates mean values of each variables during one or more cylcles
# it does not calculate crank angle resolved parameters
import sys

sys.path.append("../")


# 插值函数
def xi(theta, theta0, theta1):
    if theta < theta0 or theta > theta1:
        raise Exception("Wrong error!!")
    from math import cos, pi
    return 0.5 * (1 - cos(pi * (theta - theta0) / (theta1 - theta0)))


class IdealCylcle:
    """
    缸压重构的过程，
    self.CylGeo:气缸几何参数
    self.ValveTiming:进排气门正时

    self.mix:缸内的混合气体

    self.T2:燃烧始点的温度
    self.p2:燃烧始点的压力
    """

    def __init__(self, CylGeo, ValveTiming=None, p0=1.e5, T0=20 + 273.15, pe0=None):
        """
        初始化气缸几何参数和气阀参数
        :param CylGeo: 气缸几何参数
        :param ValveTiming: 进排气正时参数，如果没有则为默认值
        :param p0:环境压力，Pa
        :param T0:环境温度，K
        :param pe0:排气压力，正常情况下等于进气压力，考虑到潜艇等有被压存在
        """
        from ArrayTable import ArrayTable
        from Valve import ValveDesign
        self.CylGeo = CylGeo

        if ValveTiming is None:
            self.ValveTiming = ValveDesign()
        else:
            self.ValveTiming = ValveTiming

        self.p0 = p0
        self.T0 = T0

        if pe0 is None:
            self.pe0 = p0
        else:
            self.pe0 = pe0

        # 记录压缩过程
        self.CompressData = ArrayTable(5, 0)
        self.CompressData.setTableHeader(["Crank angle", "V", "p", "T", "m"])
        self.CompressData.setTableUnit(["°CA", "m^3", "Pa", "K", "kg"])

        # 记录膨胀过程
        self.ExpanseData = ArrayTable(5, 0)
        self.ExpanseData.setTableHeader(["Crank angle", "V", "p", "T", "m"])
        self.ExpanseData.setTableUnit(["°CA", "m^3", "Pa", "K", "kg"])

        # 重构的缸压
        self.Rpressure = ArrayTable(2, 0)

        # 换气过程
        self.GasExchange = ArrayTable(2, 0)

        # 最终的缸压
        self.data = ArrayTable(2, 0)
        self.data.setTableHeader(["Crank angle", "Cylinder pressure"])
        self.data.setTableUnit(["°CA", "Pa"])

    # 由充量系数和过量空气系数计算燃油消耗率，首先假设一个Ps
    def compress(self, pim=2.5e5, Tim=300, Tr=800, xr=0.0, kc=1.3, phic=1):
        """
        计算压缩过程，即从进气门关闭到燃烧开始时缸内的过程，需要估算排气温度
        :param pim:进气管压力，Pa
        :param Tim: 进气管温度,K
        :param Tr: 排气温度，K
        :param xr: 废气系数
        :param kc: 压缩多变指数
        :param phic: 充量系数
        :return:
        """
        from GasProperty import DieselMixture, Rg
        self.kc = kc
        ivc = self.ValveTiming.IVC
        self.pim = pivc = pim

        Tim = 313. + 5. / 6. * (Tim - 273.15)  # 新鲜充量温度估算经验公式
        self.Tim = Tivc = Tim * (1 - xr) + xr * Tr

        # 单缸的新鲜充量质量
        mivc = phic * pivc * self.CylGeo.V(180) / Rg() / Tivc
        self.mix = DieselMixture()

        # 由进气量、残余废气系数和循环喷油量初始化缸内气体，先随意给定
        self.mix.init_With_Mc_r(mivc, xr, mivc / 14.3 / 1.1)
        print("Intake air mass {} mg".format(mivc * 1.e6))
        print("Tivc=", Tivc)

        Vivc = self.CylGeo.V(ivc)

        # 压缩过程线
        from Valve import mod
        for i in range(ivc, ivc + 720):
            V = self.CylGeo.V(i)
            T = Tivc * pow(Vivc / V, kc - 1)
            p = pivc * pow(Vivc / V, kc)
            m = p * V / self.mix.Rg_gas(0) / T
            self.CompressData.append([mod(i), V, p, T, m])
        self.CompressData.doQuickSort(0)

        # for i in range(self.CompressData.row):
        #     self.CompressData.table[0].data[i] = mod(self.CompressData.table[0].data[i])

    # 预混燃烧过程和扩散燃烧过程
    def Burn(self, Rp, SOC=-5, alpha=1.5, Hu=42700e3, L0=14.3):

        assert 0<=Rp<1,"premixed burn fraction value error!!"
        self.Hu = Hu
        self.Rp = Rp
        self.L0 = L0
        self.alpha = alpha
        self.SOC = SOC

        # 循环喷油量
        self.mix.gf = self.mix.M_air(0) / L0 / alpha
        print("Injected fuel per cylcle per cylinder{}mg".format(self.mix.gf * 1.e6))

        self.T2 = self.CompressData.linearInterpolate(SOC, 3)
        self.p2 = self.CompressData.linearInterpolate(SOC, 2)

        # 计算预混燃烧过程
        def fun(x, T):
            return (Hu - self.mix.u(x, T)) / (alpha * L0 + x) / self.mix.cv(x, T)

        x = 0
        T = self.T2
        step = Rp / 100.
        while x < Rp:
            T += fun(x, T) * step
            x += step
        print("Temperature after premixed burn {}".format(T))

        # 计算扩散燃烧过程
        Ttemp = T
        self.p3 = self.p2 * T / self.T2
        print("Pressure after premixed burn {}".format(self.p3))

        def fun2(x, T):
            return (Hu - self.mix.h(x, T)) / (self.mix.cp(x, T) * (self.alpha * self.L0 + x))

        step2 = (1 - self.Rp) / 100.
        while x < 1:
            T += fun2(x, T) * step2
            x += step2
        print("Temperature after disfusion burn {}".format(T))

        # 计算燃烧终点的曲轴转角
        self.V3 = self.CylGeo.V(0) * T / Ttemp
        self.EOC = self.CylGeo.getFi(self.V3)
        print("End of combustion {}".format(self.EOC))
        self.T3 = T
        # self.p3=self.mix.M_total(1)*self.mix.Rg_gas(1)*self.T3/self.V3
        # print("Pressure after burned {}".format(self.p3))

        return T

    def Expense(self, ke=1.33):
        # 根据燃烧终点的状态去计算膨胀过程缸内温度和压力的变化
        self.ke = ke
        from numpy import arange
        from Valve import mod
        for i in arange(self.SOC, self.SOC + 360):
            V = self.CylGeo.V(i)
            p = self.p3 * pow(self.V3 / V, ke)
            T = self.T3 * pow(self.V3 / V, ke - 1)
            self.ExpanseData.append([mod(i), V, p, T, self.mix.M_total(1)])
        self.ExpanseData.doQuickSort(0)

    def pressureReconstruct(self, m=1):
        """
        重构缸压
        :param m: 形状因子
        :return:
        """
        from Cylinder import WibeFunction
        from numpy import arange

        # 重构缸压
        for i in arange(self.SOC, self.EOC):
            temp = WibeFunction(i, self.SOC, self.EOC - self.SOC, m=m)
            p = (1 - temp) * self.CompressData.linearInterpolate(i, 2) + temp * self.ExpanseData.linearInterpolate(i, 2)
            self.Rpressure.append([i, p])

        from numpy import arange
        for i in arange(self.SOC, self.EOC):
            self.data.append([i, self.Rpressure.linearInterpolate(i, 1)])

    def pit(self,etaTK=0.56, pik=None):
        """
        计算涡前压力
        :param pik: 压气机压比,默认为进气压力除以环境压力
        :param etaTK: 涡轮增压器效率
        :return:
        """
        if pik is None:pik=self.pim/self.p0
        from Compressor import piT
        from GasProperty import k_exhaust_gu,k_Justi
        # 计算排气门打开时的压力和温度
        pevo = self.ExpanseData.linearInterpolate(self.ValveTiming.EVO, 2)
        Tevo = self.ExpanseData.linearInterpolate(self.ValveTiming.EVO, 3)
        # pevo = self.ExpanseData.linearInterpolate(180, 2)
        # Tevo = self.ExpanseData.linearInterpolate(180, 3)
        print("Temperature at EVO {}".format(Tevo))

        # 等熵膨胀到涡前压力
        def fun(x):
            k = k_exhaust_gu(500)

            def Ttfun(kk):
                Tt0 = Tevo * pow(x * self.pe0 / pevo, (kk - 1) / kk)
                return k_Justi(Tt0, self.alpha) - kk

            h = 0.01
            while abs(Ttfun(k)) > 1.e-5:
                k -= Ttfun(k) / ((Ttfun(k + h) - Ttfun(k - h)) / 2. / h)
            # print(k)
            Tt = Tevo * pow(x * self.pe0 / pevo, (k - 1) / k)
            return piT(pik, etaTK, Tt, self.T0, self.alpha * 14.3) - x
        # self.pe0

        # 等容膨胀到涡前压力
        # def fun2(x):
        #     Tt = self.pe0 * x / pevo * Tevo
        #     return piT(pik, etaTK, Tt, self.T0, self.alpha * 14.3) - x

        x = 2.0
        h = 0.01
        while abs(fun(x)) > 1.e-5:
            x -= fun(x) / ((fun(x + h) - fun(x - h)) / 2. / h)

        # while abs(fun2(x)) > 1.e-5:
        #     x -= fun2(x) / ((fun2(x + h) - fun2(x - h)) / 2. / h)

        print("Pressure ratio of turbine {}".format(x))
        print("pressure before turbine {} bar".format(x * self.pe0 / 1.e5))
        print("Temperature before turbine {} K".format(Tevo * pow(x * self.pe0 / pevo, (1.33 - 1) / 1.33)))
        self.pem = x * self.pe0
        self.Tt=Tevo * pow(x * self.pe0 / pevo, (1.33 - 1) / 1.33)
        return self.pem

    def gasExchange(self):
        pem = self.pem
        from numpy import arange
        int = self.ValveTiming.int = (self.ValveTiming.EVC + self.ValveTiming.IVC) / 2.
        for i in arange(self.ValveTiming.EVC, int):
            self.data.append([i, self.pim])
        for i in arange(int, self.ValveTiming.IVC):
            temp = xi(i, int, self.ValveTiming.IVC)
            p = self.pim * (1 - temp) + self.CompressData.linearInterpolate(i, 2) * temp
            self.data.append([i, p])

        exh = self.ValveTiming.exh = (self.ValveTiming.EVO + self.ValveTiming.IVO) / 2.
        for i in arange(self.ValveTiming.EVO, exh):
            temp = xi(i, self.ValveTiming.EVO, exh)
            p = self.ExpanseData.linearInterpolate(i, 2) * (1 - temp) + pem * temp
            self.data.append([i, p])
        for i in arange(exh, self.ValveTiming.IVO):
            self.data.append([i, pem])

        from Valve import mod
        for i in arange(self.ValveTiming.IVO, self.ValveTiming.EVC + 720):
            temp = xi(i, self.ValveTiming.IVO, self.ValveTiming.EVC + 720)
            p = pem * (1 - temp) + self.pim * temp
            self.data.append([mod(i), p])

        # from numpy import arange
        # for i in arange(self.ValveTiming.EVC, self.ValveTiming.IVC):
        #     self.data.append([i, self.GasExchange.linearInterpolate(i, 1)])

    def plot(self):
        self.data.doQuickSort(0)

        import matplotlib.pyplot as plt
        fig, ax = plt.subplots(1, figsize=(10, 5))
        ax.plot(self.data.table[0].data, self.data.table[1].data)
        plt.xlabel(self.data.table[0].ColName + "(" + self.data.table[0].ColUnit + ")")
        plt.ylabel(self.data.table[1].ColName + "(" + self.data.table[1].ColUnit + ")")

        ax.scatter(self.EOC, self.data.linearInterpolate(self.EOC, 1))
        ax.annotate('EOC %.3g $^\circ$CA' % self.EOC,
                    xy=(self.EOC, self.data.linearInterpolate(self.EOC, 1)), xycoords='data',
                    xytext=(0, 10), textcoords='offset points',
                    arrowprops=dict(arrowstyle="->"))

        index = self.data.findMaxValueIndex(1)
        maxpreCA = self.data.table[0].data[index]

        ax.scatter(maxpreCA, self.data.linearInterpolate(maxpreCA, 1))
        ax.annotate('maxium pressure %.5g bar' % (self.data.linearInterpolate(maxpreCA, 1) / 1.e5),
                    xy=(maxpreCA, self.data.linearInterpolate(maxpreCA, 1)), xycoords='data',
                    xytext=(0, 10), textcoords='offset points',
                    arrowprops=dict(arrowstyle="->"))

        ax.scatter(self.ValveTiming.IVC, self.data.linearInterpolate(self.ValveTiming.IVC, 1))
        ax.annotate('IVC %.3g $^\circ$CA' % self.ValveTiming.IVC,
                    xy=(self.ValveTiming.IVC, self.data.linearInterpolate(self.ValveTiming.IVC, 1)), xycoords='data',
                    xytext=(0, 10), textcoords='offset points',
                    arrowprops=dict(arrowstyle="->"))

        ax.scatter(self.ValveTiming.EVO, self.data.linearInterpolate(self.ValveTiming.EVO, 1))
        ax.annotate('EVO %.3g $^\circ$CA' % self.ValveTiming.EVO,
                    xy=(self.ValveTiming.EVO, self.data.linearInterpolate(self.ValveTiming.EVO, 1)), xycoords='data',
                    xytext=(0, 10), textcoords='offset points',
                    arrowprops=dict(arrowstyle="->"))

        ax.scatter(self.ValveTiming.EVC, self.data.linearInterpolate(self.ValveTiming.EVC, 1))
        ax.annotate('EVC %.3g $^\circ$CA' % self.ValveTiming.EVC,
                    xy=(self.ValveTiming.EVC, self.data.linearInterpolate(self.ValveTiming.EVC, 1)), xycoords='data',
                    xytext=(0, 10), textcoords='offset points',
                    arrowprops=dict(arrowstyle="->"))

        ax.scatter(self.ValveTiming.IVO, self.data.linearInterpolate(self.ValveTiming.IVO, 1))
        ax.annotate('IVO %.3g $^\circ$CA' % self.ValveTiming.IVO,
                    xy=(self.ValveTiming.IVO, self.data.linearInterpolate(self.ValveTiming.IVO, 1)), xycoords='data',
                    xytext=(0, 10), textcoords='offset points',
                    arrowprops=dict(arrowstyle="->"))

        # ax.scatter(self.ValveTiming.int, self.data.linearInterpolate(self.ValveTiming.int, 1))
        ax.annotate('$p_{im}$ %.4g bar' % (self.pim / 1.e5),
                    xy=(self.ValveTiming.int, self.data.linearInterpolate(self.ValveTiming.int, 1)), xycoords='data',
                    xytext=(-28, 40), textcoords='offset points',
                    arrowprops=dict(arrowstyle="->"))
        ax.annotate('$p_{em}$ %.4g bar\n$T_t$ %.4gK' % (self.pem / 1.e5,self.Tt),
                    xy=(self.ValveTiming.exh, self.data.linearInterpolate(self.ValveTiming.exh, 1)), xycoords='data',
                    xytext=(-0, 40), textcoords='offset points',
                    arrowprops=dict(arrowstyle="->"))

        ax.scatter(self.SOC, self.data.linearInterpolate(self.SOC, 1))
        ax.annotate('SOC %.3g $^\circ$CA' % self.SOC,
                    xy=(self.SOC, self.data.linearInterpolate(self.SOC, 1)), xycoords='data',
                    xytext=(-80, 10), textcoords='offset points',
                    arrowprops=dict(arrowstyle="->"))
        #
        # ax.scatter(self.EVC, table.linearInterpolate(self.EVC, 1))
        # ax.annotate('EVC %.3g $^\circ$CA' % self.EVC,
        #             xy=(self.EVC, table.linearInterpolate(self.EVC, 1)), xycoords='data',
        #             xytext=(0, 10), textcoords='offset points',
        #             arrowprops=dict(arrowstyle="->"))
        plt.xticks([-360, -180, 0, 180, 360], ["-360\nTDC", "-180\nBDC", "0\nTDCF", "180\nBDC", "360\nTDC"])

        # ax.axhline(y=0, color='r', linestyle="-.")
        ax.axvline(x=0, color='g', linestyle=":")
        ax.axvline(x=180, color='g', linestyle=":")
        ax.axvline(x=-180, color='g', linestyle=":")
        ax.axvline(x=360, color='g', linestyle=":")
        ax.axvline(x=-360, color='g', linestyle=":")

        plt.tight_layout()
        plt.show()

    def analyze(self, speed=1500, plot=False):
        # 将压缩过程的缸压记录到重构的曲线中
        from numpy import arange
        for i in arange(self.ValveTiming.IVC, self.SOC):
            self.data.append([i, self.CompressData.linearInterpolate(i, 2)])

        # 将燃烧结束到排气门打开这一过程记录
        from numpy import arange
        for i in arange(self.EOC, self.ValveTiming.EVO):
            self.data.append([i, self.ExpanseData.linearInterpolate(i, 2)])

        from ArrayTable import ArrayTable
        from Cylinder import MeEff
        tcylce = 30 * 4 / speed

        PV = ArrayTable(3, 0)
        self.data.doQuickSort(0)
        for i in range(self.data.row):
            PV.append(
                [self.data.table[0].data[i], self.CylGeo.V(self.data.table[0].data[i]), self.data.table[1].data[i]])

        # PV.plot(2, 1)

        # 单缸指示功
        IndecatedWork = PV.integrate(2, _colx=1)

        # 单缸指示功率
        IndicatedPower = IndecatedWork/ tcylce

        # 指示平均有效压力
        IMEP = IndecatedWork / self.CylGeo.displacedVolume()

        # 指示热效率
        etait = IndecatedWork / (self.mix.gf * self.Hu)
        print("thermal efficiency {}".format(etait))

        # 迭代求解平均有效压力
        # \[IMEP*{\eta _m}\left( {{C_m},BMEP} \right) = BMEP\]
        def fun(BMEPpara):
            return BMEPpara - IMEP * MeEff(self.CylGeo.bore, self.CylGeo.stroke * speed / 30, BMEPpara)

        BMEP = IMEP
        h = IMEP / 1.e2
        while abs(fun(BMEP) > 1.e-5):
            BMEP -= fun(BMEP) / ((fun(BMEP + h) - fun(BMEP - h)) / 2 / h)
        print("BMEP={} bar".format(BMEP / 1.e5))

        etam = MeEff(self.CylGeo.bore, self.CylGeo.stroke * speed / 30, BMEP)
        print("Mechanical efficiency {}".format(etam))

        # 计算有效功率,W
        BrakePower = BMEP * self.CylGeo.totalDisplacedVolume() / tcylce
        print("Brake power {}kW".format(BrakePower/1.e3))

        # 计算油耗率,kg/h
        B = self.mix.gf*self.CylGeo.num_of_cylinders / tcylce * 3600

        # 计算燃油消耗率，g/(kW*h)，与进气质量和过量空气系数相关
        BSFC = B * 1.e3 / (BrakePower / 1.e3)
        print("BSFC={}g/(kW*h)".format(BSFC))

        print("maximum pressure {}bar".format(max(self.data.table[1].data)/1.e5))

        if plot:
            self.plot()

        return BMEP,BSFC,max(self.data.table[1].data)

    def loop(self):
        # 输出平均有效压力，油耗，
        # 气缸几何参数，充量系数，过量空气系数，CA0,燃烧形状因子
        # 调参，排气温度，进气压力，进气温度,涡轮增压器效率
        pass


if __name__ == "__main__":
    from Cylinder import CylinderGeometry

    Geo = CylinderGeometry("WP7")
    # Geo.plotVolume()
    # I = IdealCylcle(Geo)
    # I.compress(pim=1.5e5, phic=0.95,kc=1.37)
    # I.Burn(Rp=0.026, alpha=2)
    # I.Expense(ke=1.346)
    # I.pressureReconstruct(m=2)
    # I.pit(etaTK=0.41)
    # I.gasExchange()
    # I.analyze(speed=2100, plot=True)
    #
    # ps=2e5;SOC=-5;alpha=2;
    Rp=0.1
    m=2
    pim=2.5e5
    # phic=0.9

    speed=2100

    BMEPGoal=17.59e5

    # m是影响CA50的关键参数，m越小燃烧越粗暴，通常柴油机转速越低，m值越大，平均有效压力越高，m值越大，一般中、低速柴油机m，一般在0.5-1.0之间————内燃机工作过程仿真技术

    # 


    def func(phic,etaTK,m,SOC):
        I = IdealCylcle(Geo)
        I.compress(pim=pim,phic=phic)
        I.Burn(Rp=Rp,alpha=2,SOC=SOC)
        I.Expense()
        I.pressureReconstruct(m=m)
        I.pit(etaTK=etaTK)
        I.gasExchange()
        BMEP,BSFC,maxpre=I.analyze(speed=speed,plot=False)

        # return -I.analyze(speed=2100,plot=False)
        result=abs(BMEP-BMEPGoal)
        # print(result)
        return result


    from sko.GA import GA
    from sko.DE import DE
    # ga=GA(func=func,n_dim=2,size_pop=20,max_iter=20,lb=[0.8,0.4],ub=[1.1,0.7],precision=0.001)

    de=DE(func=func,n_dim=4,size_pop=10,max_iter=5,lb=[0.8,0.4,0.5,-8],ub=[1.1,0.7,1,-1])
    temp=de.run()[0]
    print(temp)
    # print(ga.run())
    # I.plot()

    import pandas as pd
    import matplotlib.pyplot as plt

    Y_history = pd.DataFrame(ga.all_history_Y)
    fig, ax = plt.subplots(2, 1)
    ax[0].plot(Y_history.index, Y_history.values, '.', color='red')
    Y_history.min(axis=1).cummin().plot(kind='line')
    plt.show()


    I = IdealCylcle(Geo)
    I.compress(pim=pim, phic=temp[0])
    I.Burn(Rp=Rp, alpha=2)
    I.Expense()
    I.pressureReconstruct(m=m)
    I.pit(etaTK=temp[1])
    I.gasExchange()
    I.analyze(speed=2100, plot=True)
    # I.plot()
